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Abstract: A comprehensive model for explaining scintillation yield in liquid xenon is intro- 
duced. We unify various definitions of work function which abound in the literature and incorporate 
all available data on electron recoil scintillation yield. This results in a better understanding of elec- 
tron recoil, and facilitates an improved description of nuclear recoil. An incident gamma energy 
range of 0(1 keV) to 0(1 MeV) and electric fields between and 0(10 kV/cm) are incorporated 
into this heuristic model. We show results from a Geant4 implementation, but because the model 
has a few free parameters, implementation in any simulation package should be simple. We use 
a quasi-empirical approach, with an objective of improving detector calibrations and performance 
verification. The model will aid in the design and optimization of future detectors. This model is 
also easy to extend to other noble elements. In this paper we lay the foundation for an exhaustive 
simulation code which we call NEST (Noble Element Simulation Technique). 
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1. Introduction and Motivations 

Liquid noble elements (LNE) have been established as an attractive detection medium in experi- 
ments searching for dark matter or neutrino-less double-beta decay: a large number of existing or 
planned experiments employ liquid xenon, ai^gon, or neon []I], 0]. The direct dark matter search ex- 
periments looking for WIMPs (Weakly Interactive Massive Particles) aim to maximize their ability 
to discriminate between electron and nuclear recoils (ER and NR respectively) within their tai^get 
volume. While the advantages of LNEs are abundant, one drawback is the non-linear energy depen- 
dence of the scintillation yield per unit energy deposited in the medium, for both ER and NR ^. 
Further, this phenomenon is influenced by the magnitude of the applied electric field [@]. In this 
article, we have addressed these issues for ER in liquid xenon by developing a general model for 
scintillation and ionization yields. We have implemented this model in a Geant4 |^ simulation and 
have demonstrated its ability to reproduce a wide variety of measurements, both for the field-free 
case and for high applied field. We have also extended our model to include NR scintillation yield 
in the zero-field case. 

An understanding of light yield is especially important in the low-energy region of interest 
(~5-50 keV) which contains the majority of the spectrum for NR from WIMPs in xenon ([Q] and 
ref. therein). The energy dependence of light yield is non-monotonic at low energies, stemming 
from effects related to dE /dx, incident particle type, and electric field. This poses a challenge for 
simulations that employ constant scintillation yields. This model provides a comprehensive Monte 
Carlo similar to that of Tawara et al. [Q] for Nal(Tl). Although this paper focuses solely on liquid 
xenon, in Section [7| we point to simple steps needed to include other noble elements. 
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Figure 1. General schematic for the ionization, scintillation, and loss processes occurring after energy 
deposition by any recoiling species in a noble element. Ionization leads to electrons either recombining 
or escaping (as described in the text). Recombination forms excimers which contribute more scintillation, 
essentially identical to that generated by direct electronic excitation (collectively, SI). Escape electrons can 
be harnessed with an electric field and their number determined with a charge amplifier coupled to an anode, 
or, in a two-phase detector (see Section 0), these electrons can produce more scintillation in a higher-electric- 
field gas phase (S2). 



2. Physical and Mathematical Framework 

We begin by reviewing the basic principles of scintillation physics. A particle depositing energy 
creates excitation and ionization, as well as heat. Initial excitation contributes to the production of 
light called the primary scintillation, or SI, through an intricate process involving excited molec- 
ular states (unnecessary to fully simulate to reproduce the light yield) [0-0] ■ The partitioning of 
deposited energy into excitation and ionization depends on incident particle type and possibly on 
electric field JTHI], but previous work indicates that the partition is independent of incident particle 
energy QTHI , [TTl ]. This entire process is summarized in Figure [I]. Some energy is lost to heat in the 
form of secondary nuclear recoils JT^ , |T3|] or sub-excitation electrons, which fail to lead to the gen- 
eration of SI light [0]. For ER, the former is negligible and the effect of the latter can be absorbed 
into a higher value for the work function (below). 

Ionization electrons either recombine and lead to additional (SI) scintillation or escape. Two- 
phase LNE direct dark matter detectors [@] drift such electrons with an electric field and extract 
them into a gas volume where they drift at higher velocity and generate more scintillation through 
more excitation and recombination in the gas. This is termed secondary, or S2, scintillation. The 
ratio of S2 to SI light differs for ER and NR ([Q, |7|] and ref. therein), and hence, it is the most 
critical discriminating parameter in dark matter searches. The first step in fully understanding this 
phenomenon is to have a comprehensive description of ER behavior in xenon. 

The distribution of the ionization electron population between recombination and escape con- 
stitutes the crux of the energy dependence of the S 1 yield. The direct excitation contribution to S 1 , 
which depends on the small ratio of excitons to ions, is consequently small (IJT^l and ref. therein). 
Since an ionization electron lost to recombination cannot escape and one which escapes cannot 
recombine, SI and S2 light yields are anti-correlated. Furthermore, scintillation yield, and hence 
recombination probability, is known to be a function of the energy loss per unit of path length 
dE/dx [§] and of the apphed electric field 
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For the energy partition we use a simplified Platzman equation [Q, |7|] 

Edep=N,,We, + NiWi=Ni{aW,, + Wi), a=NejNi (2.1) 

where E^ep is energy deposited by a particle in a single interaction (not its total energy deposition), 
Nex is the number of excitons created per deposition, Ni the number of ions, Wex the work function 
(required energy) for exciting atoms, and Wj the ionization work function. One sums the numbers 
of excitons and ions stemming from single interactions (recoils) to calculate the total yield of one 
incident particle. The ratio of excitons to ions, often labeled as a, may differ for NR versus ER, 
and may change with field, but it is likely not a function of energy JTHI]. The theoretical value 
of this dimensionless ratio a is a constant 0.06 for liquid xenon jn] , |T3p , used successfully to fit 
experimental data QT^. Some measurements of a in liquid xenon have been as high as 0.20, though 
ai^guably consistent with 0.06 within experimental uncertainties 0, [T7| ]. 

Ignoring heat, including the kinetic energy of ionization electrons (which can be absorbed 
into Wi), and assuming near- 100% efficiency for excited or recombined electrons to lead to SI 
(well-established experimentally QTHI]), the numbers of quanta produced aie 

^ Edep 

aWex + Wi 

Nex = otNi (2.3) 

which in turn lead to Nph photons and Ne ionization electrons: 

Nph=Nex + rNi and (2.4) 
Ne=Ni{\-r) (2.5) 

where r is the recombination probability. For long particle tracks we use [Q, [TE| ] 

+C, C=l-A/B, (2.6) 



dx 



derived from Birks' Law []TP|]. Doke et al. [ ]TE| ] derived a scintillation yield for LNEs from Birks' 
Law, while we adopt the same approach but instead derive the recombination probability underlying 
the scintillation yield by separating excitation from the recombination probability. The first term 
of the equation represents the recombination which occurs when a wandering ionization electron 
is captured by an ion other than its parent (termed "volume" recombination). The second term (C) 
represents "geminate" recombination, also known as Onsager recombination JTEI , PH ]. It is defined 
as ionization electrons recombining with their parent ions, which is assigned a fixed probability for 
all dE/dx. A low a implies that most S 1 light is a result of the recombination processes. 

We plot Equation in Figure with the parameter values that we will show in Section to 
reproduce the light yield of gamma rays, the first particles we consider in this work. The constraint 
is imposed so that as dE/dx approaches infinity, the recombination probability approaches unity. ^ 
Even in the absence of an electric field, some ionization electrons do not recombine on timescales 
of interest to experiments (<1 ms); they temporarily escape the Coulomb fields of the ionized xenon 
atoms and are called "escape" electrons []TE|]. For instance, ~10% of ionization electrons from a 
122 keV gamma in xenon will not be observed to produce excitons UTDI , PTj P2| ]. 



'Linear energy transfer LET, defined as the quotient of dE/dx and density, is often used in place of dE/dx. 
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Figure 2. The electron recombination probability that is implemented in our Monte Carlo simulation for 
zero electric field. It is a function of Linear Energy Transfer (LET), which is the quotient of dE/dx and 
density. We base this probability on Doke et al. [|TE|], treating A, B, and C as free parameters and using 
experimental data on light yield to constrain them. The parameters compare well with past values used such 
as C = 0.55 [^. The parameters are also constrained by the condition A /B + C = 1, which ensures complete 
recombination at infinite LET, or dE jdx. 

For short particle tracks the Thomas-Imel box model is better suited []T3|], where the recombi- 
nation probability is instead QTDI , [T3t ^ 

where a' is a constant dependent on ionization electron and hole mobilities and the dielectric 
constant (beyond the scope of this work), v is the mean ionization electron velocity, and t? is a 
length scale defining ionization density volume. We can equate a' /{a^v) with e^/{aekT) [|l0|], 
where e is the electron charge, e is the dielectric constant, k is Boltzmann's constant, and T is 
the electron temperature. Using e = 1.96, T = 165 K [P3|], and a = 4.6 |J.m [Q], we estimate a 
value of 0.14 for zero electric field. We assume here that most electrons will thermalize prior to 
recombination (and for T use the 1 atm boiling point, near which most detectors operate [@]). This 
is false at high electric field, which increases electron energy 

We define "short tracks" to be those that are shorter than the mean ionization electron-ion 
thermalization distance, 4.6 |J.m in liquid xenon [P^. If an incident gamma or electron produces 
secondary gammas or scatters multiple times, then there are multiple interaction sites. At each site, 
tracks are separately summed to determine whether to apply Thomas-Imel or Doke/Birks formal- 
ism. Each is based on the same recombination physics, albeit with different descriptions governed 
by short- or long-track limits [P5|]. The Thomas-Imel understanding relies on recombination in a 
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Figure 3. (left) A histogram of photon yields obtained in NEST in response to mono-energetic 122 keV 
gammas (^^Co). The width of the response curve demonstrates the ability of our model to generate vari- 
ation in photon yield, even when no standard means of scintillation fluctuation is included, such as that 
described by the Fano factor or detector effects. Even at this minimal level, our simulation produces 
non-Gaussianity in the low-energy tail, (right) Markedly similar real data 0]. Compare with the asymmetry 
of our simulated peak. The reason for observing spread in final light yield, sans Fano factor and detector 
effects, lies in the dependence of the recombination probability on a detailed track history and stochastic 
variations in dE /dx values. Because of ~10% photon detection efficiency, our x-axis (photons) differs by a 
factor of ten (right axis in photoelectrons). Also, the low-energy peak, in the data on the right, comes from 
partial energy deposition, which is not possible in our simulated infinite xenon volume. 



"box" (cubic) geometry [|T^, ^ while Doke constructs his formulation utilizing the cross-section 
of a long column of electron-ion pairs QTEH. We seamlessly integrate these two models into one 
unified picture of scintillation yield, which is continuous versus energy and consistent with data. 

Our model determines recombination for individual recoils in a Geant4 simulation by ex- 
plicitly tracking energies and recoil ranges for secondary electrons produced by an incident gamma. 
Doing so preserves macroscopic characteristics of the energy-dependent light yield (Figure ^ while 
remaining stochastic. As demonstrated in Figure ^ Monte Carlo fluctuations in tracks result in a 
variation in the number of scintillation photons ultimately produced, even before fluctuation in the 
ratio of excitons to ions, or their sum (Fano factor), is taken into consideration.^ 

For simplicity, our model defines a mean work function for production of either excitons or 
electron-ion pairs, such that Edep = W{Nex This is equally as effective as defining two distinct 
work functions for successfully explaining experimental results [Q, [T^ 0. Different values for the 
two processes are difficult to extract experimentally, and using a mean value causes no loss of 
generality [ [TOj , ^5] ]. Because W is also equivalent to (aW^ + Wi) / { \ + a) based on Equations ^71] 
to P3| , a higher Wgx is canceled out by a lower W,- or vice versa, leaving the total numbers of both 
quanta unchanged. Similarly, a higher a is counteracted in the above equation by a lower Wgx, 
making an exact experimental knowledge of a less relevant. 

In general, dE/dx increases with decreasing energy for electrons less than 1 MeV p5|]. Hence, 

^The latter effect is well-known to be small, sub-Poissonian in fact [pi liol UM. 
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Figure 4. (top) This figure is re-printed from [^^. It shows the anomalous decrease in scintillation per unit 
energy with increasing dE/dx (decreasing gamma energy) observed in liquid xenon for the lowest energy 
gammas studied, (bottom) A demonstration of a similar \\ig\\-dE / dx scintillation yield decrease for electron 
recoils in Nal(Tl). Here the decrease is caused by a well-known mechanism involving the kinematics of 
interactions within the crystal. This figure is re-printed from [ PIU . Note that the x-axis of the top figure is 
energy, while for the bottom it is dE /dx. 



recombination probability (and scintillation yield) should increase with decreasing energy, as long 
as the Doke model applies 0, |TE|]. However, an interesting phenomenon occurs in a low-energy 
region where the scintillation yield instead decreases (Figure For gamma rays below 0(10 keV), 
the slope of the light yield curve changes sign as it does for Nal(Tl) in a different energy 
regime [PPj, PTD, albeit due to a depletion of activation sites in the crystal [^TJl- 

One would expect recombination to be enhanced at high dE/dx due to greater ionization den- 
sity, but it has been shown that at low energies (high dE/dx) recombination becomes independent 
of ionization density QKll , [T3| , PPQ . This is accommodated in our model because low-energy parti- 
cles are strictly in the Thomas-Imel regime, where the recombination probability depends on the 
energy, via the number of ions, and not on dE/dx. In other words, for particle tracks smaller than 
the thermalization distance of ionization electrons, the total track length is no longer relevant. Dahl 
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was the first to apply a modified Thomas-Imel box model to successfully explain the odd turnover 
in Ught yield at low energy, which changes position in energy for different electric fields HTDI ] . At 
higher energies, the Thomas-Imel box model had already been successfully utilized to fit data on 
the global scale, such as total charge yield (total number of ionization electrons extracted) versus 
electric field, glossing over microphysics JT^. In Dahl's approach, the model is applied to the scale 
of the lowest energy (sub-keV and eV) individual electron recoils using PENELOPE [P^. 

In our model, we apply the Thomas-Imel model to the electron-recoil energy scale 0(1 keV) 
accessible with Geant4, making the assumption that it can be adapted for zero electric field. Man- 
alaysay et al. HT^ demonstrate that the Thomas-Imel parameterization, even though undefined at 
zero field, can have a finite limit and produce realistic results. Moreover, Figure ^ in this work 
demonstrates the utility of this assumption by successfully reproducing light yields. 



3. Zero-Field Scintillation Yield vs. Energy 



The parameters of Equation 2^ and a'/ {4a^v) of 2/7 (henceforth ^ /Nj) were treated as free in 
order to match our compilation of existing data on zero-field gamma ray scintillation. A concerted 
effort was made to include all known experimental results in our compilation and any possible 
omitted works are unknown to the authors rather than intentionally excluded. Figure ^ compares 
the results from simulations for this work to the available empirical data. 

A traditional analysis with errors could not be performed because many papers do not quote 
error bars in either their figures or their text (e.g., [030), or state that their errors are negligible [P^. 
Furthermore, there are several instances of more than one data point for a given gamma energy, with 
one or more points lacking errors, and yields are most often quoted as relative to measurements 
made at 122, 511, or 662 keV. The latter is due to experimental issues, such as difficulties in 
knowing PMT quantum efficiency at low temperatures and UV wavelengths and/or inadequate 
knowledge of the reflectivity of the walls of the detector. Thus, assumptions had to be made in order 
to translate all results into absolute yields. We used 511 and 122 keV as benchmarks, translating 
uses of 662 keV or other energies by experiments as benchmarks as needed. For the low-energy 
data of Obodovskii and Ospanov [^9|], the 59.5 keV Am line is used. 

The process of translating the data into absolute yield increases the en^or by an amount which 
can not accurately be estimated due to the varying age and quality of the data, and in some rare 
cases, due to variation in choice of benchmark energy, which has a profound effect. For instance, if 
the results of Barabanov et al. [ 031 ] ai^e re-scaled assuming a benchmark of 122 keV in Figure ^ then 
59.5 and 122 keV are not outliers, but all other points become too low. Choices like this needed to 
be made to bring the majority of a large quantity of old and new data into rough agreement. Overall, 
the ultimate consistency among different detectors is remarkable. Precedence for such comparison 
of different data sets is exemplified by Doke et al. [Q], who were forced to assume a match in yield 
between one ^^Na gamma line and a high-energy electron in order to incorporate the results of 
Barabanov et al. into their Figure 4. We treat data all as equally as possible in Figure ^ without 
averaging results except when necessary, such as when we calculated a mean 122 keV yield in 
order to incorporate data which referenced it but not 511 keV (which was our primary benchmark 
because it was as common as 122 keV in experiments, but provided the additional advantage of 
forcing the majority of data from Barabanov et al. to agree with other works). 
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Figure 5. Experimental values for absolute S 1 yields for ER in xenon as a function of incident gamma energy 
compared with our Monte Carlo output, the foundation of which is the recombination probability curve in 
Figure 0. Fitted Gaussian means (see Figure 0, left) are used to report the simulated scintillation yields. Data 
are taken from [||, [T^, [T^ H 0-^7]]. When authors quote relative yields, we infer absolute numbers 
of photons using the 511, 122, or 59.5 keV line as a reference energy. All data are translated based on the 
works where absolute yields are available [0, [T^, ^ The systematic error within one detector should be 
comparable at each energy, making relative measurements meaningful. The cluster of data points at 122 keV 
have been separated for clarity. For these and other overlapping points, some error bars are thicker than 
others for differentiation. The 3-5 MeV points are from a preprint of and are reproduced here with 
permission of M. Yamashita for completeness, though are considered less reliable due to poorer resolution. 
The 9.4 and 32.1 keV de-excitations of ^^"'Ki [ ll6| ] are not comprised solely of gamma rays, but are also 
included in the figure for completeness. Two points, at 164 and 236 keV (hollow blue crosses [P^), are from 
inelastic neutron scatters, resulting in gammas of these energies together with a nuclear recoil component, 
which may possibly serve to increase the yield. The right-hand y-axis uses the definition Wpi, = Ejep/Npi,, 
as described in the text. More features and significant outliers are addressed in Section ^ where the work of 
Doke et al. which is in good agreement with our own, is also explained in depth. 
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Table 1. The parameters of our model that best reproduce the data. Note that the value of E, /Nj is comparable 
to our earlier estimate of 0. 14. The minimum mean-squared residual for the high energy data alone, lying 
primarily in the Doke/Birks regime, was calculated in two different ways: treating multiple points for the 
same energy independently, or averaging all the points at the same energy. Since some data sets are without 
stated uncertainty, the latter was a raw, not error-weighted, average and thus less reliable. The parameters 
resulting in the best reproduction are the same in either case. 



data window and other conditions 


parameters for best reproduction (this work) 


low energy points alone (<15 keV) 
high energy points treated independently 
averaging points (raw) at same energy 


^/Ni = 0.19 
A = 0.18,B = 0.42,C = 0.57 
A = 0.18,B = 0.42,C = 0.57 



We minimize the mean-squared of the residuals in order to optimize the parameters of our 
model, which is founded on Equations and P77| . We perform a series of full Monte Carlo 
simulations with different sets of parameters until finding the optimal set. Then our model can 
successfully reproduce the available data, thus predicting yield in regions where the data is sparse. 
The two clear outliers from Barabanov et al. [ 031 ] are excluded. They may be explained by non- 
uniform light collection in the detector, as claimed by Also discarded, for the purpose of our 
plot and the optimization, was one result considerably distant from all others, quoting less than 
34 photons/keV at 22 and 88 keV Q?^. Below 15 keV, the Thomas-Imel box model dominates, 
and it is in this regime where we find the optimal value for /Ni. We then use that value for the 
data as a whole, keeping it fixed. The top of Figure ^ is a plot of the percent residuals between 
various data sets and our model. Given our assumptions for relative-absolute translation, most data 
do not differ by more than ±4% from a simulation with the parameters which best reproduce the 
data. These are listed in Table |T], above. Throughout this work, we do not fit data on scintillation 
yield directly. Instead we systematically explore ansatz solutions for the recombination probability 
which are then entered into the simulation to reproduce some measured yields and predict others. 

In addition to the list in Table [1], other parameters which we employed, but that we held fixed, 
included W = 13.7 eV, a = 0.06, and 4.6 |xm as the electron-ion thermalization length scale, which 
we pragmatically selected to determine the cross-over energy where the Doke model would begin 
to dominate, because of how well it predicts the approximate location of the major change in slope 
in scintillation yield as a function of energy. The W-value is an error-weighted average from reports 
by Doke et al. [§] and Dahl 13.8 ± 0.9 and 13.7 ± 0.2 eV. At least three other values exist, 
but were not included in the average: 14.0 eV because it was reported without a clearly stated 
uncertainty pDH, 13.46 it 0.29 eV [ P7| ] due to a calibration issue expressed by Dahl [ P5| ] (but it was 
still used to estimate a 122 keV yield), and 14.7 eV |T3|, 0, superseded in the past thirty-five 
years by more measurements benefiting from technological improvements. 

The y-axis labeled "effective" work function^ in Figure ^ is for Wp/, = Edep/Nph, which de- 
pends upon incident energy, because Nph varies with energy via the recombination probability r 
(Equation This definition expedites comparison with other works ([P5|, and ref. therein) 
which do not relate all results to one W . The ratio of W to Wph is described as scintillation effi- 



This work function is also known as W' |0] or VV^, conflicting with others' use of the same nomenclature 
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ciency'^ and the relationship Wph > W must hold true . Equality implies excitation alone, with 
no ionization. The inverse of our W = 13.7 eV corresponds to 73 photons/keV. 

The clear dip in light yield at ~30-35 keV is present because of a resonance at the xenon 
K-edge for production of low-energy Auger electrons [ 051 , 0-^. They are low enough in en- 
ergy (< 10 keV) to be governed by the Thomas-Imel model. Furthermore, at or above about 29.8 
keV, the K-shell x-ray is emitted, traveling a small distance to create a second interaction site for 
which r is separately calculated, serving to enhance the effect of the scintillation decrease: tracks 
are summed separately at each site to determine which recombination treatment to use, enhancing 
the probability of being in the Thomas-Imel regime. We note that the simulated dip is too shallow, 
possibly due to an inability of Geant4 to produce electrons of low enough energy to accurately de- 
pict the necessary interactions, even with the Auger physics module activated. In Geant4, a gamma 
does not generate electrons whose energy is so low that their range would be below the minimum 
distance cut-off. In order to maintain energy conservation, an energy deposition is recorded for the 
gamma in the tracking output file. To partially mitigate this issue, we assume that an ER with that 
recorded energy occurs [P5|]. 

Thus far we have only discussed gamma rays. We now treat electi'ons as the incident particles 
and subject them to a similar simulation as was done for incident gammas. We obtain Wph = 19.4 — 
14.7 eV, for incident electron energies in the range 1 - 100 keV. For electrons below 0(100 keV) in 
energy, Seguinot et al. report Wph = 14.3 eV p5[ ] in agreement (<5%) with the lower edge of our 
predicted range. Given the energy dependence of Wpi,, Seguinot et al.'s measurement must be an 
effective value, most likely an average or a minimum. In their earlier work PSQ, they claim a value 
of 12.5 or 12.7 eV, in disagreement with their later work and our present work, but Chepel et al. [ P5| ] 
cite Seguinot et al. as 12.7 it 1.3 eV, where the upper bound on the uncertainty (14.0 eV) disagrees 
by only ~5% with our prediction (14.7 eV). We cannot, however, reconcile their claim of a low 
ionization work function (defined later in this section) for low-energy electrons of 9.76 eV ^ 
using any of the possible definitions. Further, it has been disputed EMl- 

Unlike for low-energy electrons, data on high-energy ones are plentiful [^, derived from ^°^Bi 
decay, which is a mixture of 0.976 and 1.048 MeV electrons, and a 1.064 MeV gamma, plus other 
lower-energy decay products PPQ. Authors typically quote yield as that of an approximate "1 MeV 
electron," implying that the electrons dominate; cited yields are all in agreement with each other 
within uncertainties ^. We take two typical empirical results to compare to our simulation: 
42.0 ± 8.4 photons/keV and 46.4+43 photons/keV.^ The model as it stands predicts too high a 
yield, 57.4 photons/keV for a 1 MeV electron. As an aside, for a 1 MeV gamma the simulated yield 
of our model is 59.2 photons/keV. This difference in yield is explained by the lower recombination 
probability for the ionizations from the ER of an incident 1 MeV electron compared to the lower- 
energy (higher-dE / dx) ER produced by a (Compton scattering) 1 MeV gamma ray. However, this 
difference is too small to explain the much lower yield in 1 MeV electron data. 

Our initial attempt to match the absolute yield of 1 MeV electrons using the Doke model was 
not successful. Forcing a match results in too low of an absolute yield for other well-established 
results. Setting the (relative) 122 keV point from Barabanov et al. [ P3| ] equal to one of the absolute 

'^This efficiency should not be confused with the light collection efficiency of a detector. 
^0.64 ± 0.03 relative yield with respect to VK = 13.8 eV [Q]. 
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Table 2. The predictions of our model for the S 1 photon yield of a 1 MeV electron at zero electric field, 
under various changes in parameters. Our model can reproduce the reported experimental yield to high 
precision, but only under extreme conditions, fully explained below in the text. 



conditions (parameters) 


yield (ph/keV) 


A = 0.18,C = 0.57,a = 


0.06 (unchanged model parameters) 


57.4 


A = 0.00,C = 0.57,a = 


0.06 (geminate/Onsager recombination alone) 


43.1 


A = 0.00,C = 0.57,a = 


0.20 (additional initial excitation) 


46.7 


most recent experimental data, with the smallest uncertainty 


46.4 



measurements of the 122 keV yield lowers the yield of the 1 MeV gamma, hence that of the 1 MeV 
electron as well, creating better agreement with data. Unfortunately, this makes other Barabanov 
et al. points too low, in direct conflict with much more recent measurements [Q, ^ and barely 
within the quoted eiTor of Ni et al. [^. (Ni et al. report absolute, not relative, yields, requiring 
no assumptions in interpreting their results.) The problem at hand is subtly evident in the work of 
Doke et al. (the basis for the curve in Figure ^ labeled Doke 2001). They use data from Yamashita 
et al. at the expense of ignoring Barabanov et al. [ P3| ] and determine the absolute yields of 1 
MeV electrons and gammas to be nearly identical ; using 122 or 5 11 keV as the key to translate 
relative into absolute yield makes the absolute yield too high again for a 1 MeV electron as in our 
own work. In Ni et al. [Q], one possible solution can be found: using a universal (average) dE /dx 
for interactions (explained in Section the Doke model can be made to fit the 1 MeV electron 
yield along with equal or lower-energy gammas, but we wish to avoid losing the stochastic nature 
of the simulation. Finally, even if one were to try applying an existing model other than Doke's, 
the similarity in dE / dx between a gamma and an electron at 1 MeV and the direct contradiction 
between one older set and two or three newer data sets would remain points of contention. 

An extreme solution to the challenge is to assume that non-geminate recombination ceases 
for a minimally ionizing particle. Parameter A of Equation represents such recombination, and 



is proportional to ionization density, which is in turn proportional to dE/dx [jT8|]. Perhaps low 
ionization density prevents ionization electrons from recombining with non-parent ions because 
of distance, though Mozumder claims that this is impossible [P^. The thermalization distance 
(4.6 \xm) is much greater than the reach of the Coulomb field of a xenon ion as defined by the 
Onsager radius (49 nm) [Q] , so Onsager recombination should not be dominant. For the sake of 
argument, we surmise that a thermalized ionization electron can eventually be re-attracted to its 
parent ion if no other ion is available. Alternatively, while thermalizing, an ionization electron 
could become trapped in a closed path within the Coulomb field of the parent ion if repelled by 
neighboring electron clouds []5D|]. This may help explain the slow recombination time observed by 
Doke and Masuda p2|]. However, we do not present our solution as a definitive one, given a lack 
of fundamental plausibility, but instead as a pragmatic one, allowing the creation of a simulation 
which produces reasonable results for gammas and electrons alike in a wide range of energies. A 
summary of our results for the 1 MeV electron under parameter variation is presented in Table 0. 
Changing a to correspond with the results of [Q] in addition to setting A = matches data best. A 
higher a may imply more Onsager recombination than encoded in C. More such recombination is 
hard to distinguish from a shift in the initial ratio of excitons to ions []T^. 
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After reproducing zero-field scintillation yield it is possible to make predictions of measure- 
ments other than yield. Aprile et al. report that the zero-field escape fraction for ionization electrons 
from a 662 keV gamma is ;^ = 0.22 ±0.02 [|l7p (also called jS ||Sp). From our simulation we derive 
X = 0. 19. With an a of 0.20 and the adjustment of A = explained above, we can reproduce the 
X = 0.43 reported for ^°^Bi decay [Q]. In addition to x we can derive a value for the effective 
ionization work function, defined as W^-^^^ = Ejgp/Ni. The definition of W^^^^ = W( 1 + a) is equiv- 
alent 0, ^ 0].'' Our model indicates that it should be 14.5 eV, in agreement with the 14.2 ± 0.3 eV 
reported by Obodovskii and Pokachalov (as cited by [J7|]), but not with the oft-cited 15.6 ± 0.3 eV 
([^, ^ and ref. therein). The discrepancy is easily explained: 15.6 eV corresponds with a = 0.14 
if IV = 13.7 eV, or a = 0.06 coupled with the outdated value ofW = 14.7 eV [|T5|]. One more result 
we have now is our simulated ratio of W^"-^-^ to the band gap energy in liquid xenon of 1.56, in good 
agreement with the theoretical (1.65) and measured (~1.6) values [@]. We are continuing to collect 
more data to compare with the model, making improvements to the model whenever possible. 

4. Discussion of Other Models 

In this section, we discuss the advantages of our model as compared to others that are available 
in the literature. Figure ^includes the calculation of Doke et al. This approach is based on 
extracting an electron response from a portion of the gamma data [^, ^ , applying that extracted 
response to early-generation secondary electrons and then reconstructing the gamma response. In 
contrast, our model, which uses Equations and P77| , is applied to each new electron created in 
the cascade, as low in energy as permitted by Geant4. The range of each electron is obtained from 
Geant4 output, and not assumed from a Bethe-Bloch calculation. Doing so allows for variations 
in dE/dx. Such variations contribute to the spread in scintillation yield (Figures and This 
approach can be incorporated in any simulation which outputs lists of energy depositions, and 
utilized for predicting charge yields (Section ^ and the NR quenching factor (Section Unlike 
some other previous works [Q, ^, we do not utilize a direct fit to the energy dependence of total 
scintillation yield. We thus predict, not extrapolate, the response to an energy regime that is lower 
than what has been studied experimentally. 

In earlier work, similar models have been applied to materials other than LNE, for instance, 
Nal(Tl), which is a well known scintillator. The recombination probability (r) for Nal(Tl) has the 
same functional form PTD as Equation albeit with C = 0, leaving A = B. Despite the compo- 
sitional differences between mono-atomic LNE and solid crystal scintillators (often comprised of 
multiple elements, including dopants), or organic scintillators for which Birks' Law was derived 
originally, this formula has been applied to LNE with modest success IQ, ^, [TE|]. However, its ap- 
plication involved defining a universal dE / dx for a gamma interaction, which had to be obtained 
by averaging over all interactions (ERs). Electron numbers and energies of course vary [Q], dimin- 
ishing the utility of a mean dE / dx by leading to a static recombination probability. Our treatment 
of xenon avoids usage of such an average. 

There was, however, no certainty that Equation P?^ , which works relatively well with global 
dE/dx, would also work when accounting for individual dE/dx. We explored numerous other 

^Just like the other work functions, the notation is inconsistent. This one is variably called simply IV |^ or 
Wg ^.OmW is at times denoted Wph{max) [Q], W,,, [0], or W" [|l|, g]. 
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possibilities with the same qualitative features (low probability at low dE /dx and high proba- 
bility at high dE /dx), but these were unmotivated and made for poor reproductions of the data 
{erf{x), exponential, Equation with dE /dx raised to different powers, and non-unity r at in- 
finite dE /dx.) No benefit was seen in inventing a new approach, and using C = (no geminate 
recombination) did not minimize the mean-squared error, while a non-zero C did. Numerous other 
physically-motivated models exist, Mozumder's for example [P^, but have more free parameters, 
are inapplicable at zero field, with no simple adaptation, or apply only to MeV energies [P5|]. We 
achieve agreement with the most data and the fewest parameters. There also appeared to be no 
benefit to the introduction of heat loss for ER, to explain the trend of low-energy gamma data. The 
dE/dx of a low-energy ER is compai^able to that of an NR [P^, but the Lindhard mechanism, by 
which NR lose energy to additional nuclear recoils, does not apply to low-energy gammas, which 
should exclusively produce ER QTDI , [T^ . 

Our model is ultimately an approximation because Geant4 does not accurately simulate elec- 
tron production (tertiary, quaternary generations, etc.) below 250 eV [^, unlike other software 
packages with more precise low-energy electron physics, such as PENELOPE [ PlU , 0. However, 
we demonstrate that working within the default Geant4 energy regime is sufficient for all practical 
purposes. As shown in Figure ^ (top), our work can match data in most cases to better than 5%, 
assuming an accurate conversion from relative to absolute gamma and electron yield at zero field. 
In addition, the model continues to work well for non-zero field, as shown in the next section. 

There are various extensions to the existing simulation packages which address scintillation 
directly. For example, RAT [ pTl ] is a Geant4 add-on that simulates the scintillation yield of many 
elements, with the exception of xenon. Our modeling of ionization and our focus on xenon, both 
lacking in RAT [ p^ , make our work complementary. 

5. Scintillation Quenching with Electric Field 

As electric field strength increases, ionization electrons liberated by an interaction are increasingly 
less likely to recombine with an ion. This is known as electric field scintillation quenching (|@] 
and ref. therein), not to be confused with the quenching factor for NR, which involves different 
scintillation loss mechanisms such as the Lindhard effect ^ [T^ |T3|]. For a two-phase detector, this 
results in an increase in S2 light output at the expense of S 1 . The problem is quite complex because 
recombination probability does not change uniformly across all energies []T^. Dahl demonstrates 
that a fundamental theory of recombination as a function of electric field is possible, which would 
explain the trend in recombination probability using a modified Thomas-Imel box model []T^. For 
reasons delineated below, we approach the problem semi-empirically, adapting as free parameters 
A, B, and C from Equation p^ , and /Ni from Equation p77| . We tune these on a subset of available 
data sets, and then accurately predict many others. 

An alternative would be to attempt a first-principles approach in Geant4 by creating look-up 
tables for the recombination probability r based on a Thomas-Imel application in PENELOPE. For 
three reasons, we decided against this approach. First, the introduction of tables would eliminate 
the benefits of the stochastic nature of simulations, an attribute that was a primary motivation for 
implementing our model. Even a look-up table with stochastic variations fails to capitalize on the 
details of track history captured by Geant4, thus making it less realistic. Second, the first-principles 
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Table 3. Summary of the Thomas-Imel parameter /Ni in our model, as compared to Dahl (who does not 
discuss zero field). The different columns for Dahl are, respectively, ER only, both ER and NR, and ER and 
NR where a is field-dependent for NR. Our work is for ER and a constant a. We maintain a — 0.06 at all 
fields. 





Electric Field (V/cm) 


This Work 


Dahl 2009 [ 


10 


] 




60 ±5 
522 ± 23 
876 ± 36 
1951 ± 86 
4060 ± 190 


0.19 
0.034 
0.029 
0.026 
0.023 
0.021 


0.0339(2) 0.0348(2) 0.0386(2) 
0.0335(3) 0.0354(2) 0.0342(2) 
0.0296(3) 0.0303(2) 0.0299(1) 
0.0371(5) 0.0348(3) 0.0319(2) 
0.0317(5) 0.0318(4) 0.0280(2) 



approach breaks down at zero field, where single-phase detectors operate. Third, it is not verified 
experimentally above 0(100 keV) as opposed to the approach employed by Doke et al. [Q, ^, [TE|]. 
Energies in the 0.5 - 3.0 MeV range are important for double-beta decay experiments and use of 
xenon in PET scans [p3|]. Thus, our hybrid model of Thomas-Imel plus Doke covers a large range 
in both electric field and energy. 

We began our data-matching with one data set that is the most comprehensive in terms of 
simultaneous span in energy (2-200 keV) and electric field (60-4060 V/cm) QIDI]. By starting with 
one very complete data set, we initially sidestepped the issue related to collecting data sets from 
different detectors with distinct systematic effects. First, we let <^ /Nj vary, using data below 15 keV 
to determine the best reproductions of data. Our results, in Table 0, differ from the values used by 
Dahl QTDI ] for the same data; despite the disagreement the data is reproduced well, as demonstrated 
later in Figure Agreement with data in spite of disagreement in the value of a parameter from the 
same model is likely a consequence of not applying the Thomas-Imel formula to the lowest energy 
(< 250 eV) electrons. Geant4 does not generate them, and we are forced to utilize the Thomas- 
Imel box model in an approximate fashion. (We adapted Thomas-Imel to zero field at low energy 
without an extrapolation, but with a recalculation, in Section 0). Another contributing factor is the 
fact that Dahl applies the Thomas-Imel model completely to all his data, while we use it only at 
low energy, switching to Doke/Birks for the longer tracks of particles at higher energies. The lack 
of tracking of the shortest tracks causes the Thomas-Imel model to break down in Geant4 at high 
energy, but we take advantage of the regime where it does work. After having found values for 
t, /Ni at five different fields, we fit a power law to it: t, /Ni = 0.057£"° i2 (field in V/cm). This is 
consistent with Dahl's expectation of t, /N oc E^'^-^ [|l0|]. 

The next step was to find the Doke parameters as a function of electric field magnitude. We 
began this study by looking at the well-established electric field dependence of the 122 keV gamma. 
The result for 122 keV has been repeated several times ^HX with duplicate and triplicate 

points often taken at the same electric field to verify repeatability We proceeded to check the 
variation of the Doke parameters achieved studying only this one energy against other energies and 
to adjust our fit for the electric field dependence of the parameters accordingly, though continuing 
to give 122 keV results the greatest weight. 
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Figure 6. (left) Choice of parameters (A, B, and C from Equation [2.61 ) and their behavior as a function of 
electric field. There are two distinct cases. The red curve was determined by matching available gamma 
data, giving the most weight to the 122 keV gamma. In the second case, the blue curve is for relativistic 
electrons, matched to 1 MeV electron data. We propagated our extreme assumption of A = 0, and then varied 
C. Decreasing values of A and B imply decreasing volume recombination; a decrease in C corresponds to 
a decrease in Onsager recombination (see text). We use a power law fit to interpolate A and C between 
discrete fields. The final outcome is a decreasing recombination with increasing field, as expected. The case 
of zero field is a discontinuity where A = 0. 18, B = 0.42 (both zero for high-energy electrons), and C = 0.57 
(Figure and Section 0). (right) Variations in recombination probability, as a function of applied electric 
field, used in our simulation. This set of curves is based only on using C = and the variation in A (= B) 
as depicted at left. The high-energy electron case is not plotted because it is the trivial case of constant 
recombination probability with respect to LET. 



There must be a transition between the Thomas-Imel regime and the Doke regime. At zero 
field we use the ionization electron thermalization length scale as the cross-over distance, and this 
immediately makes the transition smooth. As electric field increases a need arises to change this 
cross-over distance in order to avoid creating a discontinuity. We offer a power law fit as a function 
of electric field as a first effort without optimization: the distance in |j.m is 69 /VE, where the 
electric field magnitude is in units of V/cm. 

Monotonic variation emerged naturally in our model (Figure ^ for A and B, and a negligible 
value for C, so we set C = 0. These parameters have physical meaning as discussed earlier. We in- 
terpret a near-zero value for C as demonstrating that volume recombination dominates over Onsager 
recombination. This corresponds exactly with the findings of Dahl, who reports that Onsager re- 
combination should be negligible QTOI]. With an external electric field overwhelming the local fields 
of ions, this is expected. In the case of zero field, the opposite appears true (C = 0.57), but this 
is empirically supported in argon QTEI], with available evidence being ambiguous for xenon ^. 
This is one approach, but we were unable to explain as much data as easily with a fully physically- 
motivated model for electric field, the Jaffe model, r = I — l/{l + K/E), where K is the so-called 
recombination coefficient and E is the electric field. However, the Jaffe model came close to 
matching significant amounts of data as modified by Aprile et al. to take into account delta elec- 
trons [|7|, 0. Our approach to the electric field dependence has the advantage of adapting the 
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Figure 7. (top row) Comparison of our model with experimental data on relative scintillation yield as a 
function of electric field. The relative yields of 1.0 are unique to each energy and are obtained by normalizing 
to zero-field values. In several plots our curves lie on top of each other if too similar in energy. Our model 
is valid to higher fields as well, but we restrict these plots to select ranges which are rich in data. Empirical 
data are from [|,|Tg,|T7|,0gg,g|,^B|31], labeled by the first author in each case, 
(bottom row) Comparison of our model with data on relative charge yield as a function of electric field, with 
data taken from [Q, |I^ |T7|, ^ ^ Q ^3|]. In this case, the relative yields of 1.0 are again unique to 
each energy, but correspond to the infinite-field ionization electron yield. This normalization assumes that 
at infinite field all possible charge is extracted. Thus the relative yield can equivalently be thought of as the 
escape fraction for ionization electrons separated from ions. "All possible charge" is defined as the number 
of ions or electrons generated ab initio, assuming a = 0.06 (0.20 for 1 MeV electrons). Experimental 
results quoted are appropriately re-scaled if assuming a different value of a either explicitly or implicitly 
(Wl^J~f = 15.6 eV may imply a different a for example) [ |l7| , as are ones plotting results in terms of 
the fraction of total quanta instead of the fraction of electrons [PHQ. Two contradictory 570 keV gamma 
charge yield measurements are shown. Interestingly, the data set (Aprile 1991 p9|]) that is inconsistent with 
our model matches well with the 554 and 481 keV electrons, which are associated decay products of the 
570 keV gamma from ^"^Bi decay. The green curve in the bottom right plot has been generated using this 
interpretation. Underestimated charge yield for the 368 keV electron (blue line compared to pink crosses) 
may be explained by a detector effect: ionization electron multiplication in charge collectors [W^- 
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Figure 8. The logarithm of the ratio of the number of ionization electrons (the basis of S2) to the number 
of S 1 scintillation photons, generated by the ER stemming from a high-energy gamma ray (356 keV in both 
simulation and data), at an electric field of 876 V/cm. Actual data, from the Xed liquid xenon detector, is 
depicted as a color map (red to yellow scale), reproduced from Figure 5.1 of C.E. Dahl's Ph.D. Thesis [[Till ]. 
The output of our simulation model is superposed as a simple density plot in blue. The green curve is a raw 
average of our simulated points in 1 keV bins. Vertical marks at ^30 keV are likely related to the xenon K- 
edge, or are binning artifacts. Because we do not simulate the Xed geometry, we do not study and compare 
the density of our points to the density of the data as encoded by the coloration. To approximate the Xed 
result we look at the energy deposition in a finite volume of xenon of comparable dimensions to Xed but 
without walls. The scatter in the simulation is a result of recombination fluctuations discussed in the text: 
particle history variation and dEjdx fluctuations, which swamp Fano factor and a effects. The experimental 
data have additional scatter, either due to a low-yield tail in S2 or a high-yield tail in SI. The former arises 
from detector effects, wherein ionization electrons can be absorbed by walls JTHI]. The upper dashed curve 
represents the centroid of the ER band as shown in the original figure [[T^], while the lower one is the NR 
band, beyond the scope of this work. 



existing formulae for zero field instead of introducing more models. 

Our total compilation of comparisons is available in Figure 0. The data points are shown 
without any error bars for either the yield or the electric field magnitude, for clarity. (There are 
few cases where they are even available.) A significant amount of recent data is predicted accu- 
rately, having begun by reproducing 122 keV points and the low-energy regime in Dahl's data. 
Discrepancies exist, but were inevitable, due to contradictory results being present in such a large 
compilation (likely due to non-obvious and unreported systematic errors), and due to the fact that 
we never directly fit any of these experimental yield-dependence curves, even at 122 keV. Instead, 
we reproduce or predict data with simulation by testing differing underlying recombination proba- 
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bility curves. With a few parameters that describe the recombination process for individual recoils 
and are functions of electric field, we can effectively explain global aspects such as total light and 
charge yields. These yields combine individual yields from all recoils spawned by one incident 
particle. This is a clear break from the past work of others who directly fit the total yields with 
models (for example, \^ ^ ^). 

We omit works ([]55|] for example) which define "total" light and charge yields for partial 
energy deposition "peaks," as these would require proper modeling of the detector in question in 
order to reproduce well. We also exclude older works with measurements that cannot be reconciled 
with more recent ones, including [p7|], and Voronova et al. as cited in [p^. More recent experiments 
have benefited from technological advancements, leading to more accurate measurements. Recent 
9.4 keV data is poorly predicted by the Thomas-Imel model. Curiously, agreements can be achieved 
by approximating the ^^"jjq. de-excitation as only a gamma ray and applying the Doke model 
despite the low energy. Without such adjustments these data contradict Dahl [ ]TD| ] in this energy 
range. The contradiction may arise from inaccurate decomposition into constituent electrons and 
gammas in Geant4. 

In spite of these few faults, we have successfully reproduced many measurements based on 
various detectors. Figure || shows our prediction for \ogiQ{Ne /Nph) as superposed on the data from 
Dahl IJT^ . While the low-energy regime is a reproduction based on these data, the extension above 
~15 keV is the prediction from our model. Here we have compared simulation with data at only 
one of Dahl's electric field choices, as an example. Comparable agreement is seen at other values 
of electric field as well. 

6. Nuclear Recoil Quenching Factor 

We have applied the framework motivated in Section 0, and validated in Section 0, to zero-field 
NR scintillation yield. NR scintillation efficiency has been a subject of debate, with several con- 
tradictory experimental results and theoretical underpinnings [ pEj |59| ]. This is an issue of great 
significance for dark matter detectors. Dark matter particles are expected to produce nuclear, not 
electron, recoils in most theories 

We offer a unique approach here (applied earlier only by Dahl []TD|], but at non-zero field) that 
uses pai^ameters from a model motivated by gamma-induced ER and applies them to NR directly. 
Tracks of NR below 0(100 keV) are always within the Thomas-Imel regime []T^. We again use 
the single free parameter from the Thomas-Imel model (§ /Ni) that successfully explained the low- 
energy yield of ER. We make only one modification to the mathematical framework for NR, such 
that Equation ^TT] is modified [|T5 |T3|] as follows: 

Ed,pL{Edep) = {Ne. +Ni)W = (1 + a)NiW (6.1) 

where L is the Lindhard factor (or, effective Lindhard factor for modification to the Lindhard theory 
made to better explain data). It should not be confused with ^eff, which represents the ratio 
of scintillation light produced by an NR to that produced by a 122 keV gamma, at zero electric 
field. Both are functions of energy, and the latter is the traditional way to report results on NR. 
Though not equivalent to L, knowledge of ^eff enables an empirical determination of L, which 
represents the fraction of energy available for excitation and ionization. Unlike ER, NR will lose 
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most of its energy (~80%) to heat, which in this case implies interactions with other nuclei instead 
of electrons. Nuclear interactions with electrons excite or ionize atoms and thus produce usable 
signals for a detector [JT}, |T3p . The fractional factor L is the reason why two energy scales are often 
defined: keVg^ and keV,„. (or keV, ) |@, [T^l- The former is for ER and the latter for NR. In this 
paper we use only "keV" as the total energy of a recoil as reported by Geant4. Without knowing 
L or ^eff one can not a priori define the NR energy scale in terms of the "electron equivalent" 
energy in "keVgg." 

We use the same value of /Ni (= 0.19) as for ER to evaluate the model against empirical 
data. A choice remains in the application of the Lindhard factor with respect to proposed coitcc- 
tions. We elect to investigate two scenarios in order to treat this intricate subject with due diligence: 
pure Lindhard, and the Hitachi correction to the Lindhard factor [1T3|]. Lindhard theory is more ap- 
propriate for solid crystal scintillators, and for xenon nuclei Hitachi has performed a first-principles 
recalculation QTB]]. We thus closely follow the methodology of Sorensen and Dahl, but we differ by 
starting with an electron-recoil-motivated Thomas-Imel parameter at zero electric field instead of 
one found from fits to NR data []T3|]. Therefore, we offer a novel approach. 

Our two predictions in Figure ^ are consistent with recent measurements by Plante et al. [ ]5D| ] 
and Manzur et al. |]8p , but in direct conflict with the older claims of Sorensen et al. and Aprile et 
al. [pEQ. Because NR yields change very weakly with field [^, we did not initially investigate 
field dependence. Dahl's study suggests an explanation for this weak dependence: a ~ 1 for NR at 
non-zero field and perhaps is field-dependent, not a fixed 0.06 as in ER JTHI]. (Energy dependence 
at low energy is also possible QH]]). If a significant fraction of interactions is due to excitations, 
then recombination probability is less important. For zero field we saw no need to change the a 
or t, /Ni used successfully eai^lier. To reproduce field dependence, which is of greater importance 
now as S2 is used to lower the threshold of xenon detectors []65|], we suggest employing a /Ni 
that changes differently with electric field than it does for ER, rather than attempting a change in 
a. Though this is in apparent contradiction with Dahl's work, we offer it as another means to the 
same end. The ability to distinguish between microphysics interpretations (changing a or changing 
/Ni) is likely beyond the reach of macroscopic empirical measurements that observe yields. 

7. Conclusions 

We have presented a coherent and comprehensive approach towards understanding the scintilla- 
tion and ionization processes in liquid xenon. Starting with an ansatz-driven approach, aided by 
physically motivated models, we have been able to not only explain a large majority of the world 
data, but make predictions in regions where measurements are sparse. Our model is especially 
applicable to the low recoil energy regimes of interest for direct dark matter searches, as demon- 
strated by our prediction of ^eff- NEST will be made available to the world community to be used 
as an add-on to Geant4. Changing a few parameters in NEST, not all free (W , a, A, etc.), will 
make it possible to simulate noble elements other than xenon. We have also suggested a universal 
nomenclature (W, W^v, W,-, Wp/,, W,-^^^) in order to unify different definitions in the literature, which 
would enable the community to compare and contrast results which are often quoted in varying 
and confusing fashions. Following in the footsteps of other recent works striving for a unified view 
UTIll , |T31 , P7| ] , we suggest the use of a mean work function (simply labeled as W) as the standard for 
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Figure 9. (top row) Our predictions for NR compared with all past data and predictions. For clarity, we have 
divided the data sets into discrete (left) and continuous (right), as presented by the authors. Data are gathered 
from existing compilations as well as from the primary sources [Q, Errors bands for data in the 

right plot have been omitted for further clarity. The axis labeled relative scintillation efficiency {^ef f) refers 
to a traditional normalization to ^^Co, a mixture of 122 keV (dominant) and 136 keV gammas. Note that our 
curves are not fits but predictions. We plot the same predictions at top left and top right. Surprisingly, given 
the expected inapplicability of the Lindhard theory to liquid xenon, unmodified it appears to agree the best 
with the latest measurements. It is in good agreement with Plante et al. [^] and mostly resides within 1 ff of 
Manzur et al. [@], which seems better described with the Hitachi correction. (Note; we do not include works 
that are primarily commentary on older data [ pH 0|, but do include other theoretical and phenomenological 
predictions, extrapolations, and calculations, like that of Collar [pP|].) 

(bottom) We plot the escape electron fraction for further verification of our model. Since the escape electrons 
cannot be collected at zero field, the points listed as Manzur et al. are not a direct measurement: they are an 
extrapolation based on a model Manzur et al. use to provide a good fit to their scintillation yield data, one 
which included the basic assumption of light-charge anti -correlation [^. 



reporting results and as the benchmark for quoting relative yields. We have provided a framework 
for simulating electron and nuclear recoils at different energies and applied electric fields, and have 
confirmed the results against a large sample of data from literature. 
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